﻿using System;
using System.Text;
using System.Drawing;
using System.Buffers;
using System.Collections;
using System.Collections.Generic;
using System.Runtime.InteropServices;

public static partial class NativeAOT
{
    [UnmanagedCallersOnly(EntryPoint = "treanor")]
    public static unsafe void treanor(double t, double h, int n, IntPtr y_ptr, double eps, IntPtr f_x_ya_n_da_ptr)
    {
        double* y = (double*)y_ptr.ToPointer();
        f_x_ya_n_da = Marshal.GetDelegateForFunctionPointer<delegatefunc_x_ya_n_da>(f_x_ya_n_da_ptr);

        treanor(t, h, n, y, eps);
    }

    /// <summary>
    /// 变步长Treanor方法
    /// f计算微分方程组中各方程右端函数值的函数名。
    /// </summary>
    /// <param name="t">积分起始点。</param>
    /// <param name="h">积分步长。</param>
    /// <param name="n">一阶微分方程组中方程个数，也是未知函数个数。</param>
    /// <param name="y">存放n个未知函数在起始点t处的函数值。返回n个未知函数在t+h处的函数值。</param>
    /// <param name="eps">控制精度要求。</param>
    public static unsafe void treanor(double t, double h, int n, double* y, double eps)
    {
        int i, j, m;
        double x, dt, hh, pp, s, aa, bb, dd, g, dy, dy1;
        double* w = stackalloc double[4 * n];
        double* q = stackalloc double[4 * n];
        double* r = stackalloc double[4 * n];
        double* d = stackalloc double[n];
        double* p = stackalloc double[n];
        double* u = stackalloc double[n];
        double* v = stackalloc double[n];

        hh = h;
        m = 1;
        pp = 1.0 + eps;
        x = t;
        for (j = 0; j <= n - 1; j++)
        {
            u[j] = y[j];
        }
        while (pp >= eps)
        {
            for (j = 0; j <= n - 1; j++)
            {
                v[j] = y[j];
                y[j] = u[j];
            }
            t = x;
            dt = hh / m;
            for (i = 0; i <= m - 1; i++)
            {
                for (j = 0; j <= n - 1; j++)
                {
                    w[j] = y[j];
                }
                f_x_ya_n_da(t, y, n, d);
                for (j = 0; j <= n - 1; j++)
                {
                    q[j] = d[j];
                    y[j] = w[j] + h * d[j] / 2.0;
                    w[n + j] = y[j];
                }
                s = t + h / 2.0;
                f_x_ya_n_da(s, y, n, d);
                for (j = 0; j <= n - 1; j++)
                {
                    q[n + j] = d[j];
                    y[j] = w[j] + h * d[j] / 2.0;
                    w[n + n + j] = y[j];
                }
                f_x_ya_n_da(s, y, n, d);
                for (j = 0; j <= n - 1; j++)
                {
                    q[n + n + j] = d[j];
                }
                for (j = 0; j <= n - 1; j++)
                {
                    aa = q[n + n + j] - q[n + j];
                    bb = w[n + n + j] - w[n + j];
                    if (-aa * bb * h > 0.0)
                    {
                        p[j] = -aa / bb; dd = -p[j] * h;
                        r[j] = Math.Exp(dd);
                        r[n + j] = (r[j] - 1.0) / dd;
                        r[n + n + j] = (r[n + j] - 1.0) / dd;
                        r[3 * n + j] = (r[n + n + j] - 1.0) / dd;
                    }
                    else
                    {
                        p[j] = 0.0;
                    }
                    if (p[j] <= 0.0)
                    {
                        g = q[n + n + j];
                    }
                    else
                    {
                        g = 2.0 * (q[n + n + j] - q[j]) * r[n + n + j];
                        g = g + (q[j] - q[n + j]) * r[n + j] + q[n + j];
                    }
                    w[3 * n + j] = w[j] + g * h;
                    y[j] = w[3 * n + j];
                }
                s = t + h;
                f_x_ya_n_da(s, y, n, d);
                for (j = 0; j <= n - 1; j++)
                {
                    q[3 * n + j] = d[j];
                }
                for (j = 0; j <= n - 1; j++)
                {
                    if (p[j] <= 0.0)
                    {
                        dy = q[j] + 2.0 * (q[n + j] + q[n + n + j]);
                        dy = (dy + q[n + n + n + j]) * h / 6.0;
                    }
                    else
                    {
                        dy = -3.0 * (q[j] + p[j] * w[j]) + 2.0 * (q[n + j] + p[j] * w[n + j]);
                        dy = dy + 2.0 * (q[n + n + j] + p[j] * w[n + n + j]);
                        dy = dy - (q[n + n + n + j] + p[j] * w[n + n + n + j]);
                        dy = dy * r[n + n + j] + q[j] * r[n + j];
                        dy1 = q[j] - q[n + j] - q[n + n + j] + q[n + n + n + j];
                        dy1 = dy1 + (w[j] - w[n + j] - w[n + n + j] + w[n + n + n + j]) * p[j];
                        dy = (dy + 4.0 * dy1 * r[n + n + n + j]) * h;
                    }
                    y[j] = w[j] + dy;
                }
                t = t + dt;
            }
            pp = 0.0;
            for (j = 0; j <= n - 1; j++)
            {
                dd = Math.Abs(y[j] - v[j]);
                if (dd > pp) pp = dd;
            }
            h = h / 2.0;
            m = m + m;
        }
        return;
    }

    /*
    // 变步长Treanor方法例
      int main()
      { 
          int i,j;
          void tnrf(double,double [],int,double []);
          double t,h,eps,y[3];
          y[0]=1.0; y[1]=0.0; y[2]=-1.0;
          t=0.0; h=0.001; eps=0.0000001;
          cout <<"t = " <<t;
          for (i=0; i<=2; i++)
              cout <<"  y(" <<i <<") = " <<setw(10) <<y[i];
          cout <<endl;
          for (j=1; j<=10; j++)
          { 
              treanor(t,h,3,y,eps,tnrf);
              t=t+h;
              cout <<"t = " <<t;
              for (i=0; i<=2; i++)
                  cout <<"  y(" <<i <<") = " <<setw(10) <<y[i];
              cout <<endl;
          }
          return 0;
      }
    // 计算微分方程组中各方程右端函数值	
      void tnrf(double t, double y[], int n, double d[])
      { 
          t=t; n=n;
          d[0]=-21.0*y[0]+19.0*y[1]-20.0*y[2];
          d[1]=19.0*y[0]-21.0*y[1]+20.0*y[2];
          d[2]=40.0*y[0]-40.0*y[1]-40.0*y[2];
          return;
      }
    */
}

